function [f, entrant, lambda_e, equity_mat, debt_mat, lambda_mat, coupon_mat, optimal_firm, optimal_coupon, C] = equilibrium(bParams,oParams,setts,upMat)
% Calculates the equilibrium level of creative destruction f

% Variables:
% bParams,oParams,setts	see readme
% upMat		Updating matrix equity value

% Set parameters
    tax = oParams.pi;
    H = bParams.H * (1 - tax) * (1- setts.EntrySubsidy);
	f_prec = setts.f_prec;
	f_max = setts.f_max;
	f_min = setts.f_min;
    f_max_iterations = setts.f_max_iterations;
	nodebt = setts.nodebt;
	debt_static = setts.debt_static;

%% Solve no debt case
	setts.nodebt = 1;
	setts.debt_static = 0;

% Set initial value error f
	f_error = f_prec + 1;

% Fixed point iteration using binomial search
	while f_error > f_prec

% Set f as middle point interval
		f = (f_max + f_min)/2;

% Calculate the entrants equity value
		[entrant, lambda_e, equity_mat, debt_mat, lambda_mat, coupon_mat, optimal_firm, optimal_coupon, C] = entrant_value(bParams,oParams,setts,upMat,f);
        
% Update the bounds for f
		if entrant < H
			f_max = f;
		else
			f_min = f;
        end

% Update the error
		f_error = abs(entrant-H);
	end

% Display result
	disp(['Solved no debt case and f is ',num2str(f)])

%% Solve debt case

% Reset settings
	setts.nodebt = nodebt;
	setts.debt_static = debt_static;

	if nodebt==0
% Set initial value error f and other parameters
		f_error = f_prec + 1;
		f_max = f+0.01;
		f_min = f;
		f_dum = 0;

% Set intial value f_history and entrant_history
        f_history = [];
        entrant_history = [];
    
% Fixed point iteration using binomial search
		while f_error > f_prec

% Save time
			t0 = clock;

% Set f as middle point interval
			f = (f_max + f_min)/2;
            
% Save value f and entrant
            f_history = [f_history; f];
            entrant_history = [entrant_history; entrant];
            
% Update recovery value
            setts.recovery_value = optimal_firm;

% Calculate the entrants equity value
			[entrant, lambda_e, equity_mat, debt_mat, lambda_mat, coupon_mat, optimal_firm, optimal_coupon, C] = entrant_value(bParams,oParams,setts,upMat,f);
            
% Update the bounds for f (where f_max is increased until for f_max entrant < H)
			if entrant < H
				f_max = f;
				f_dum = 1;
			else
				if f_dum == 0
					f_max = f_max+0.025;
					f_min = f;
				else
					f_min = f;
				end
            end

% Update the error
			f_error = abs(entrant-H)/H;
			disp([' Entrant minus H = ',num2str(f_error)]);
			disp([' Time for a single solution equity value = ',num2str(etime(clock,t0)/60),' minutes']);
            
% Regression to find f if it doesn't converge
            if length(f_history)==f_max_iterations
                X = linspace(min(f_history),max(f_history),10000);
                Herror = (H-interp1(f_history,entrant_history,X)).^2;
                [~,index] = min(Herror);
                f_max = X(index);
                f_min = f_max;
                scatter(f_history,entrant_history,"*")
            end

% Set error to zero after calculating the results for the regression estimate f
            if length(f_history)>f_max_iterations
                f_error = 0;
                disp(' Algorithm did not converge so I used linear interpolation to find a best estimate');
            end
                
		end

% Display result
		disp(['Solved debt case and f is ',num2str(f)])

	end

end
